Wannier functions analysis of the nonlinear Schrodinger equation with a periodic 

potential 



O 

o 



On 



O 

in 
o 

(N 
O 



c 

o 
o 



G. L. Alfimov*, P. G. Kevrekidis**, V. V. Konotop* and M. Salerno* 
* F. V. Lukin's Institute of Physical Problems, Zelenograd, Moscow, 103460, Russia. 
** Department of Mathematics and Statistics, 
University of Massachusetts, Amherst, MA 01003-4515, USA. 
t Departmento de Fisica and Centro de Fisica da Materia Condensada, Universidade de Lisboa, 
Complexo Interdisciplinar, Av. Prof. Gama Pinto 2, Lisboa 1649-003, Portugal 
* Dipartimento di Fisica "E.R. Caianiello" , Universitd di Salerno, 1-84081 Baronissi (SA), Italy, ana 
Istituto Nazionale di Fisica della Materia (INFM), Unitd di Salerno, Italy. 

(February 1, 2008) 

In the present Letter we use the Wannier function basis to construct lattice approximations of the 
nonlinear Schrodinger equation with a periodic potential. We show that the nonlinear Schrodinger 
equation with a periodic potential is equivalent to a vector lattice with long-range interactions. 
For the case-example of the cosine potential we study the validity of the so-called tight-binding 
approximation i.e., the approximation when nearest neighbor interactions are dominant. The results 
are relevant to Bose- Einstein condensate theory as well as to other physical systems like, for example, 
electromagnetic wave propagation in nonlinear photonic crystals. 
PACS numbers: 42.65.-k, 42.50. Ar,42.81.Dp 



Interplay between nonlinearity and periodicity is the 
focus of numerous recent studies in different branches of 
modern physics. The theory of Bose-Einstein conden- 
sates (BEC) within the framework of the mean field ap- 
proximation [ffl] is one of them. Recent interest in the 
effects of periodicity in BEC's has been stimulated by 
a series of remarkable experiments realized with BEC's 
placed in a potential created by a laser field Q (the so- 
called optical lattice). Nonlinearity and periodicity have 
been observed to introduce fundamental changes in the 
properties of the system. On the one hand periodicity 
modifies the spectrum of the underlying linear system 
resulting in the potential of existence of new coherent 
structures, which could not exist in a homogeneous non- 
linear system. On the other hand, nonlinearity renders 
accumulation and transmission of energy possible in "lin- 
early" forbidden frequency domains; this, in turn, results 
in field localization. This situation is fairly general and 
can be found in other applications, such as the theory 
of electromagnetic wave propagation in periodic media 
(so-called photonic crystals) ^ . 

The study of nonlinear evolution equations with pe- 
riodic coefficients is a challenging and interdisciplinary 
problem. This problem cannot be solved exactly in the 
general case and thus gives rise to various approximate 
approaches. One of them, borrowed from the theory of 
solid state is the reduction of a continuous evolution 
problem to a lattice problem (i.e., reduction of a partial 
differential equation to a differential-difference one). It 
turns out that the relation between the properties of pe- 
riodic and discrete problems is indeed rather deep (for 
a recent discussion of the relevant connections see e.g., 
H and references therein). Following the solid state ter- 
minology here we will refer to a discrete approximation 
when only nearest neighbor interactions are taken into 



account as a tight-binding model. This model has re- 
cently been employed in the description of BEC in an 
optical lattice || . One of the advantages of the lattice 
approach is that it allows one to obtain strongly local- 
ized configurations, the so-called intrinsic localized modes 
(ILM) (also called breathers) 0, in a rather simple way. 
These entities correspond to gap solitons of the origi- 
nal continuum model. In the above mentioned works a 
formal analysis has been provided, using a basis of func- 
tions strongly localized about the minima of the periodic 
potential. This basis, however, has not been presented 
explicitly and even its existence has not been established. 

In this work we propose to use Wannier-function (WF) 
j|,H as a complete set of functions localized near the min- 
ima of the potential, to reduce the evolution of a nonlin- 
ear partial differential equation with periodic coefficients 
to a nonlinear lattice. WF has recently been used both in 
connection with BEC in optical lattices P| and in connec- 
tion to gap solitons in nonlinear photonic crystals JT3] . In 
our case this approach leads to a vector set of lattice equa- 
tions. These lattice equations exactly correspond to the 
original continuum problem and the scalar tight-binding 
approximation can be deduced from them under some 
specific conditions. Checking these conditions one can 
analyse the applicability of the tight-binding model. In 
particular, we argue that although the ILM's reported in 
H do exist, their dynamics and stability must be stud- 
ied within the framework of a more general vector-lattice 
equation. 

Being interested in BEC applications we base our 
analysis on the ubiquitous example of the nonlinear 
Schrodinger equation 
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where a = ±1 and V(x) is a periodic potential V(x+L) 
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V(x) [fL4[ . Consider the eigenvalue problem associated 
with 



dx 2 



+ V{x)tpk,a = E a (k)tpk,, 



(2) 



where ipk, a has Bloch (Floquet) functions (BF's) tpk, a — 
e lhx Uk, a (x), with Uk, a (x) periodic with period L, and a 
is an index which labels energy bands E a (k). As is well 
known, E a (k + ^) = E a (k); thus one can represent 
the energy as a Fourier series 



n,a — ^ na 



< a (3) 



where an asterisk stands for complex conjugation and 
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The BF's constitute an orthogonal basis. However, for 
our purposes it is more convenient to use the WF's in- 
stead of the BF's. We recall that the WF centered around 
the position nL (n is an integer) and corresponding to the 
band a is defined as 

r~L f r / L 

w a (x-nL) = \ — ip k , a {x)e- lnkL dk. (5) 
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Conversely, 



^ w n>a (x)e 



in kL 
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Similarly to BF's, they form a complete orthonormal 
(with respect to both n and a) set of functions, which, 
by properly choosing the phase of the BF's in (^), can 
be made real and exponentially decaying at infinity 
In what follows we assume that this choice is made: 
w na( x ) = w n, a {x). Due to completeness of WF's, any 
solution of (□) can be expressed in the form 



ip(x,t) = 2_,Cn, a {t)w n>a {x) 



(7) 



which after substitution in (|l]) gives 
dr 



dt 
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where 



W nn in2 n3 _ I j 
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are overlapping matrix elements. Since WF's are real, 
W^2 2 a 3 at * s symmetric with respect to all permuta- 
tions within the groups of indices (a, a\, a^, as) and 



(n,ni,n2,ri3). Eq. (||) can be viewed as a vec- 
tor discrete nonlinear Schrodinger (DNLS) equation for 
c n =col(c„i, c„2, ■■•) with long-range interactions. In its 
general form, Eq. (|J) is not solvable; however it allows 
reductions to simpler lattices in a number of important 
special cases. Below we list some of them. 

(i) If the coefficients of the Fourier series (S) decay 
rapidly and |cDi ja | 3> |a) n ,a|, n > 1 one can neglect long- 
range interaction terms in the linear part of Eq. (||) tak- 
ing into account nearest neighbors only. 

(ii) Since w nia (x) is localized and centered around 
x = nL, one can assume that in some cases among all 
the coefficients W^a 2 al those with n — n\ = n 2 = n% 
are dominant and other terms can be neglected. Then, 
taking into account points (i), (ii) one arrives at the equa- 
tion 
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which degenerates into the tight-binding model 



dr 



dt 
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if one restricts consideration to the band a only. Note 
that within the single band approximation, Eq. ( |Tl| ) can 
be generalized by including next nearest neighbor over- 
lapping terms from Eq. ^, thus leading to the mixing 
of on-site and intra-sitc nonlincaritics of the same type 
as in the model introduced in [ fU| |. It should also be 
mentioned that the coefficients W™™a 2 a 3 m Eq. (0) are 
independent of n. 

(iii) In the general case, however, single band de- 
scriptions can become inadequate (see below) due to 
resonant interband interactions induced by nonlinearity 
(this is quite different from linear solid state physics 
where interband transitions are usually induced by ex- 
ternal forces). In this case Eq. (|lfj| ) can be further 
simplified by supposing that the periodic potential de- 
pends on some parameter e: V{x) = V e (x), such that 
cji q, = £)i ta (e) — O(e) when e — * 0. After the trans- 
formation c nia (t) = exp{iu)o ta t}c nta (t) one arrives at 
the equation for c n , a with explicit dependence on t in 
the nonlinear terms in the form of oscillating exponents 
exp[i(a) , Q + U) , ai - w , Q2 - &o, a3 )t\. Let also c„, Q (0) 
be small enough. Then on the timescale 1/e these ex- 
ponents are rapidly oscillating unless a = 0:2,^1 = ct3 
or a = a 3 ,ai = a 2 Then, denoting W™™ ai = W aai 
(the coefficients W aai do not depend on n and describes 
interband interactions), and using time averaging tech- 
niques Q, one can reduce the lattice equation ( |l0| ) to 
the form 
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1 ~ — — a [C n -i a + C n +i a ) + 



dt 



+ <y Warn I 
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This is a vector DNLS equation with coupling between 
bands of the cross phase modulation type To inves- 
tigate ILM solutions in the Wannier representation we 
can restrict to the scalar case described by Eq. ([ll]) for 
which construction of ILM's is well-established 0. ILM's 
with multiple components of c n _ a populated can also be 
constructed (see below). 

Several comments about the above assumptions are 
in order. Firstly, the latter imply that the procedure 
of reduction of the NLS with periodic coefficients to a 
lattice is a multistep process, and thus different lattices 
will appear for different regions of the parameters. Sec- 
ondly, for the reduction to be consistent, the parameters 
of the problem must provide us with a small parame- 
ter. Thus the largest of the quantities u) n ,a/^>i,a (n > 1) 
and W£}%2*/W2ZZ£ ia9 (n 3 £ n) will define this small 
parameter of the problem. This, in particular, means 
that simplification of the lattice equation, and hence the 
reasoning for the reduction to a lattice model, are (po- 
tentially) not always available for all parametric regimes, 
and must be verified for each model. 

In the present Letter we study the validity of the above 
assumptions for Eq. (Q) with the potential V(x) = 
Acos(2x) (which corresponds to the typical experimental 
setting for BEC in optical lattices Q|). In this case Eq. 
(||) is the Mathieu equation. Table I shows the coeffi- 
cients <l> n ,a for the three lowest energy bands for A = —1 
and A = -15. 





A = -l 


A = -15 


n 


a = 1 


a = 2 


a = 3 


a = 1 


a = 2 


a = 3 





0.1305 


2.4657 


6.3604 


-9.7862 


0.0421 


8.4659 


1 


-0.1428 


0.5426 


-1.0067 


-0.0005 


0.0151 


-0.1798 


2 


0.0204 


0.0784 


0.0529 


0.0000 


0.0001 


0.0091 


3 


-0.0048 


0.0481 


-0.1107 


-0.0000 


0.0000 


-0.0008 


4 


0.0014 


0.0225 


0.0140 


0.0000 


0.0000 


0.0001 



TABLE I. The first five Fourier coefficients i2i n a of the low- 
est energy bands for two values of the amplitude potential. 

In general, one can conclude that the greater \A\ is, the 
better the linear part obeys the nearest neighbour ap- 
proximation, which is intuitively expected since the prob- 
ability of tunclling between neighbor potential wells de- 
creases with the amplitude of the potential. At the same 
time, if A is fixed the coefficients Lj n ,a, n = 0, ±1, .. de- 
cay faster for lower bands a. The results illustrate, that 
the nearest neighbor approximation works for both po- 
tential amplitudes, while the averaging resulting in ( O ) 
is applicable for A — —15 but not for A = — 1. The rea- 
son is that in the latter case the frequencies of oscillating 
exponents in (iii) are of order of uii >a . 



Moving to assumption (ii), let us introduce the fol- 
lowing notation. We denote by A r " m (A) the number of 
coefficients W^ll 2 2 ^ 3 , N < n, Oj < m, i,j = 1,2,3 (the 
coefficients with permuted indices are regarded as differ- 
ent) such that |W°2iS5S|l > A - As ii; is clear A P la y s 
the role of the small parameter of the second condition, 
and N™ m (A) gives the number of sites/zones necessary 
to take into account for maintaining the given accuracy. 
In the cases of the amplitudes A = — 1 and A = —15 
we have obtained that iV^O.l) = ATfi(O.Ol) = 1 for 
n = 1, 5. For AT£ m (0.1) and AT™ m (G.dl), see Table II. 



A 


n 


7V™ 2 (0.1) 


7V™ 2 (0.01) 


^3(0.1) 


A^yo.oi) 


-1 





4 


4 


7 


13 




1 


4 


48 


7 


219 




2 


4 


54 


7 


249 




3,4 


4 


60 


7 


303 




5 


4 


60 


7 


339 


-15 





4 


4 


13 


14 




1-5 


4 


4 


13 


26 



TABLE II. The values AT" m (A). 



Table III presents the overlapping coefficients W aai for 
two values of the amplitude of the cos-like potential. 



A 


W n 


W 22 


W33 


W 12 


Wis 


1^23 


-1 


0.375 


0.240 


0.173 


0.182 


0.152 


0.142 


-15 


0.892 


0.623 


0.473 


0.417 


0.262 


0.326 



TABLE III. Overlapping coefficients W aai ■ 

It follows from Tables II and III that, in general one 
cannot neglect the contribution of the WF of the highest 
zones. However, one can show that the model (|ll| ) can be 
successfully used to describe bright monochromatic GS 
solutions of ([l]) of the form ip(t,x) = e lut u{x). An exam- 
ple is shown in Fig. |] for (the "intermediate" between 
the above presented ones case of) A = — 5. The two pan- 
els show the cases of w = —1.5 (left panels) and oj = 1.5 
(right panels), for a = 1. The top panels show the com- 
parison of the exact solution (shown by solid line) of Eq. 
(Q) with the reconstructed profile obtained from solving 
Eq. (|2|) and using Eq. p). The relevant profiles in the 
tight binding approximation are shown by dashed line, 
while in the right panel (where the one band approxima- 
tion is less accurate), the 3-band approximation is also 
shown by dash-dotted line. The bottom panels show in a 
semilog plot the square modulus of the configurations of 
the top panels as well as, additionally, by dotted line the 
result of time evolution (for t w 50) of Eq. (|l|) with the 
tight binding approximation as the initial condition of 
the simulation. One can straightforwardly observe that 
the approximate solution "reshapes" itself into the exact 
solution (possibly shedding some very small amplitude 
radiation wakes in the process). This demonstrates that 
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the method can be used very efficiently to construct (ap- 
proximate) solutions of the original PDE, by using the 
lattice reduction in the WF representation. 



Taking into account (a)-(d), making the substitution 







FIG. 1. Comparison of the lattice reconstructed solution in 
the tight binding (dashed line) and the 3-band (dash-dotted 
line) approximation with the exact solution (solid line). The 
comparison is performed for A = —5 and u) = —1.5 (left 
panels) and A = —5, uj = 1.5 (right panels). In the bottom 
(semilog) panels additionally the result of dynamical time evo- 
lution of the tight binding approximation is shown by the dot- 
ted line. The latter can be seen to approach, as time evolves 
(the shown snapshots are for t « 50), the shape of the exact 
solution (in the left panel it can actually not be distinguished 
from it) and to match its asymptotics, possibly shedding small 
wakes of low amplitude wave radiation in the process (see e.g., 
the bottom right panel). 

Let us return now to the requirement (iii) and argue 
that choosing the small parameter as e = one can 



provide averaging of (Q) in the limit A 
we claim that 

(a) If a is fixed and A — > — oo then 



-oo. Namely, 



u ,q ~ A + (2a - l)\/ = ! + ((2a - l) 2 + l)/8 . 

If A is fixed then u>o. a tends to infinity as a grows; 

(b) If a is fixed then the value u>i iCt tends to zero faster 
than any power of 1/|j4|; at the same time if A is fixed 
then (Di iQ tends to infinity as a grows; 

(c) If a is fixed and A — > -co the Wannier functions 
can be approximated by the formula 



(2\A\)i 



TT^^2 a ~ 1 {a- 1)! 



W\Y< 



where Hk(y) are Hermite polynomials. This is a natural 
consequence of the fact that for sufficently low levels the 
potential can be well approximated by the parabolic one; 

(d) The coefficients W™™%* with different n,m,n 2 
and ri3 tend to zero as A —* — oo and at the same time 
w aaia 2 a 3 ~ K aaia2Ct3 \A\i where K a . aua2 , a3 do not de- 
pend on A and can be expressed explicitly through the 
integrals of products of Hermite polinomials [12] . 



,«(*) 



-1/8 = 



(t) and averaging over rapid 



oscillations one arrives at ( p_2[ ) with W aai — K aaiaai . 

To conclude we have shown how to derive lattice mod- 
els which approximate efficiently nonlinear partial differ- 
ential equations with periodic coefficients. This analysis 
gives the possibility to control the validity of the tight- 
binding approximation. In particular, we have shown 
that in a large region of parameter space, for the cos-like 
potential, one cannot restrict consideration to the low- 
est band. This is due to interband transitions originat- 
ing from the nonlinearity (a situation very different from 
the one known in the (linear) solid-state physics, where 
the interband transitions occur due to effect of pertur- 
bations). However, there exist parameter ranges where 
with reasonably high accuracy the atomic wave function 
(that is a bright gap soliton of the ID NLS equation) is 
approximated by a single WF. Such a state will form a 
"Wannier-soliton" that should also be experimentally ob- 
servable. It should be highlighted that the use of the WF 
basis allows one to test, extend and improve the tight- 
binding approximation, in a controllable and systematic 
fashion by accounting for higher order terms in the Wan- 
nier expansion. Moreover, there is a computational gain 
when computing with a discrete system with respect to 
the corresponding cost for a much finer mesh (needed to 
resolve the original continuous system). While this gain 
may not be overly significant in one dimension, it may 
prove quite useful in tackling higher dimensional prob- 
lems. 

It should be stressed that even though developed for a 
specific, physically relevant (to optical lattices in BEC) 
setting, the approach presented here is very general and 
directly applicable to numerous other physical problems 
including the description of solitary wave propagation 
through one-dimensional photonic crystals, [16] chemi- 
cal reactions on periodic catalytic substrates, [ 17 1 or even 
population dynamics in appropriately heterogeneous sub- 
strates. |§. 
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